
refract_time=50;
l_win=1000;
r_win=5000;
min_dur=10;

do_layer=2;

cell_cnt=0;
exp_cnt=0;
snip_cnt=0;
snipS_cnt=0;
snipE_cnt=0;
vel_snip_cnt=0;
vel_snip_c_cnt=0;
mo_snip_cnt=0;
vs_snip_cnt=0;
vs_stop_snip_cnt=0;
act_ind=[];
cell_id=[];
active_cells=[];
xcorr_delay=[];
max_act_ind=[];
ps_snips=zeros(l_win+r_win+1,120000);
ps_snipsS=zeros(l_win+r_win+1,20000);
ps_snipsE=zeros(l_win+r_win+1,20000);
if exist('mo_snips','var');
    mo_snips(:)=0;
    vs_snips(:)=0;
    vs_stop_snips(:)=0;
    saccade_snips(:)=0;
else
    mo_snips=zeros(l_win+r_win+1,150000);
    vs_snips=zeros(l_win+r_win+1,250000);
    vs_stop_snips=zeros(l_win+r_win+1,100000);
    saccade_snips=zeros(l_win+r_win+1,300000);
end

ps_snip_time=[];
vel_pb_session=[];
vel_snip=zeros(15001,5000);

pcf=[];

saccade_cnt=0;



for ind=1:length(AV_data)
    disp(num2str(ind))
    if AV_data(ind).layer==do_layer & AV_data(ind).ACode==1
        exp_cnt=exp_cnt+1;
        nbr_cells(exp_cnt)=size(AV_data(ind).raw_act,2);
        type=(AV_data(ind).velP_ind&AV_data(ind).velM_ind)+2*(AV_data(ind).velP_ind&~AV_data(ind).velM_ind)+3*(~AV_data(ind).velP_ind&AV_data(ind).velM_ind);
        %type=circshift(type,1+round((length(type)-1)*rand));
        
%         type_offsets=find(diff(type>0&type<3)==-1);
%         
%         for knd=1:length(type_offsets)
%             type(type_offsets(knd):type_offsets(knd)+100)=0;
%         end
        
        conc_onsets=intersect(find(diff(AV_data(ind).velP_ind)==1),find(diff(AV_data(ind).velM_ind)==1));
        conc_onsets=conc_onsets(conc_onsets>refract_time);
        use_onsets=logical([]);
        for knd=1:length(conc_onsets)
            use_onsets(knd)= sum(AV_data(ind).velP_ind(conc_onsets(knd)-refract_time:conc_onsets(knd)))==0 & sum(AV_data(ind).velM_ind(conc_onsets(knd)-refract_time:conc_onsets(knd)))==0;
        end
        mo_start=10*conc_onsets(use_onsets);
        
            
        [ps_raw]=find_ps_during_running(AV_data(ind).velP_ind,AV_data(ind).velM_ind);
        [ps_start]=parse_cond_starts(ps_raw,min_dur,refract_time,1);
        
        %ps_start=find(diff((AV_data(ind).PS)>0.75)==1);
        
        %[mo_start]=parse_cond_starts_vs(AV_data(ind).velP_ind,min_dur,refract_time,1);
        [vs_start]=parse_cond_starts_vs(~AV_data(ind).velP_ind&AV_data(ind).velM_ind,min_dur,refract_time,1);
        
        % this will look for pb offsets during non-running
        [vs_raw]=find_pb_during_sitting(AV_data(ind).velP_ind,AV_data(ind).velM_ind);
        [vs_stop]=parse_cond_starts_vs(vs_raw,min_dur,refract_time,2);

        interim_runnig=-1*ones(4,1);
        
        pb_onsets=[];
        if length(AV_data(ind).shutter_open)==2
            pb_onsets(1)=AV_data(ind).shutter_open(1);
            pb_onsets(2:4)=find(diff(AV_data(ind).PS(AV_data(ind).shutter_open(2):end)>0.1)>0)+AV_data(ind).shutter_open(2);
            pb_offsets=find(diff(AV_data(ind).PS(AV_data(ind).shutter_open(2)+1000:end)>0.1)<0)+AV_data(ind).shutter_open(2);
            pb_dur=pb_offsets(1)-pb_onsets(2);
            interim_runnig(3)=sum(AV_data(ind).velP_ind(round(pb_offsets(1)/10):round(pb_onsets(3)/10)));
            interim_runnig(4)=sum(AV_data(ind).velP_ind(round(pb_offsets(2)/10):round(pb_onsets(4)/10)));
        elseif length(AV_data(ind).shutter_open)==3
            pb_onsets(1)=AV_data(ind).shutter_open(1);
            pb_onsets(2:4)=find(diff(AV_data(ind).PS(AV_data(ind).shutter_open(3):end)>0.1)>0)+AV_data(ind).shutter_open(3);
            pb_offsets=find(diff(AV_data(ind).PS(AV_data(ind).shutter_open(3)+1000:end)>0.1)<0)+AV_data(ind).shutter_open(3);
            pb_dur=pb_offsets(1)-pb_onsets(2);
        elseif length(AV_data(ind).shutter_open)>=4
            pb_onsets=AV_data(ind).shutter_open+1000;
            pb_dur=pb_onsets(2)-pb_onsets(1)-2000;
        else
            pb_onsets=[];
        end
        pb_onsets=pb_onsets+1000;
        pb_dur=pb_dur-2000;
        
        if length(AV_data(ind).shutter_open)==2
            tmp_PS=round(3*AV_data(ind).PS);
            tmp_PS(1:AV_data(ind).shutter_open(2))=1;
            fb_onsets=find(diff(tmp_PS==0)==1);
            fb_offsets=find(diff(tmp_PS==0)==-1);
        else
            fb_onsets=[];
            fb_offsets=[];
        end
        
        if length(fb_onsets)>3
            disp('Warning, fb_onsets too long')
        end
        
        P_snips=[];
        M_snips=[];
        for jnd=1:length(pb_onsets)
            try
                P_snips(:,jnd)=AV_data(ind).velP(pb_onsets(jnd):pb_onsets(jnd)+pb_dur);
                M_snips(:,jnd)=AV_data(ind).velM(pb_onsets(jnd):pb_onsets(jnd)+pb_dur);
            end
        end
        tmp=[];
        for jnd=2:size(P_snips,2)
            [~,tmp(jnd-1)]=max(xcorr(P_snips(:,jnd),M_snips(:,jnd)));
        end
        tmp=tmp-length(P_snips)-1;
        PM_delay(ind)=mean(tmp);
        
        for jnd=1:length(ps_start)
            try
                if sum(AV_data(ind).velP_ind(round(ps_start(jnd)/10)-50:round(ps_start(jnd)/10)))>25
                    vel_snip_cnt=vel_snip_cnt+1;
                    vel_snip(:,vel_snip_cnt)=AV_data(ind).velP(ps_start(jnd)-5000:ps_start(jnd)+10000);
                    vel_snip_ind(:,vel_snip_cnt)=AV_data(ind).velP_ind(round(ps_start(jnd)/10)-500:round(ps_start(jnd)/10)+1000);
                    vel_pb_session(vel_snip_cnt)=find(sign(ps_start(jnd)-pb_onsets)>0,1,'last');
                end
            end
        end
        
        for jnd=1:length(ps_start)
            try
                try_again=1;
                while try_again
                    rand_start=round(max(ps_start)*rand);
                    if sum(AV_data(ind).velP_ind(round(rand_start/10)-50:round(rand_start/10)))>25
                        vel_snip_c_cnt=vel_snip_c_cnt+1;
                        vel_snip_c(:,vel_snip_c_cnt)=AV_data(ind).velP(rand_start-5000:rand_start+10000);
                        try_again=0;
                    end
                end
            end
        end
        
        saccade_times_tmp=AV_data(ind).iframe_times(AV_data(ind).saccade_frames);
        use_saccades=ones(length(saccade_times_tmp),1);
        for jnd=1:length(saccade_times_tmp)
            if min(abs(saccade_times_tmp(jnd)-pb_onsets+1000))<1000
                use_saccades(jnd)=0;
            end
        end
        saccade_times_tmp=saccade_times_tmp(logical(use_saccades));
        
        
        for knd=1:size(AV_data(ind).raw_act,2)
            if max(AV_data(ind).raw_act(:,knd))<100
                cell_cnt=cell_cnt+1;
                std_noise=sqrt(sum((AV_data(ind).raw_act(AV_data(ind).raw_act(:,knd)<1,knd)-1).^2/sum(AV_data(ind).raw_act(:,knd)<1)));
                active_cells(cell_cnt)=sum(AV_data(ind).raw_act(:,knd)>(1-3.72*std_noise))>20;
                
                [ROI_act_ms]=map_ROI_act_to_aux_data(AV_data(ind).raw_act(:,knd),AV_data(ind).frame_times,length(AV_data(ind).PS));
                max_act_ind(cell_cnt)=max(ROI_act_ms);
                sta_cnt=snip_cnt+1;
                for jnd=1:length(ps_start)
                    try
                        if min(ROI_act_ms(ps_start(jnd)-l_win:ps_start(jnd)+r_win))>0
                            snip_cnt=snip_cnt+1;
                            ps_snips(:,snip_cnt)=ROI_act_ms(ps_start(jnd)-l_win:ps_start(jnd)+r_win);
                            ps_snip_time(snip_cnt,1)=ps_start(jnd);
                            ps_snip_time(snip_cnt,2)=jnd;
                            pb_session=find(sign(ps_start(jnd)-pb_onsets)>0,1,'last');
                            ps_snip_time(snip_cnt,3)=ps_start(jnd)-pb_onsets(pb_session);
                            ps_snip_time(snip_cnt,4)=pb_session;
                            ps_snip_time(snip_cnt,5)=mean(AV_data(ind).velP(ps_start(jnd)-500:ps_start(jnd)));
                            ps_snip_time(snip_cnt,6)=interim_runnig(pb_session);
                            ps_snip_time(snip_cnt,7)=ps_start(jnd)-mo_start(find(ps_start(jnd)-mo_start>0,1,'last'));
                            ps_snip_time(snip_cnt,8)=mean(AV_data(ind).velP(ps_start(jnd):ps_start(jnd)+3000));
                            ps_snip_time(snip_cnt,9)=ind;
                            ps_snip_time(snip_cnt,10)=mean(AV_data(ind).velP(ps_start(jnd)-400:ps_start(jnd)-200));
                            ps_snip_time(snip_cnt,11)=mean(AV_data(ind).velM(ps_start(jnd)-400:ps_start(jnd)-200));
                            ps_snip_time(snip_cnt,12)=mean(AV_data(ind).velP(ps_start(jnd):ps_start(jnd)+200));
                            ps_snip_time(snip_cnt,13)=mean(AV_data(ind).velP(ps_start(jnd)+501:ps_start(jnd)+1500));
                            ps_snip_time(snip_cnt,13)=cell_cnt;
                            ps_snip_time(snip_cnt,14)=knd;
                            
                        end
                    catch
                        snip_cnt=snip_cnt-1;
                    end
                end
                end_cnt=snip_cnt;
                for jnd=1:min(5,length(ps_start))
                    try
                        if min(ROI_act_ms(ps_start(jnd)-l_win:ps_start(jnd)+r_win))>0
                            snipS_cnt=snipS_cnt+1;
                            ps_snipsS(:,snipS_cnt)=ROI_act_ms(ps_start(jnd)-l_win:ps_start(jnd)+r_win);
                        end
                    end
                end
                for jnd=max(1,length(ps_start)-5):length(ps_start)
                    try
                        if min(ROI_act_ms(ps_start(jnd)-l_win:ps_start(jnd)+r_win))>0
                            snipE_cnt=snipE_cnt+1;
                            ps_snipsE(:,snipE_cnt)=ROI_act_ms(ps_start(jnd)-l_win:ps_start(jnd)+r_win);
                        end
                    end
                end
                
                if mean(mean(ps_snips(1:1000,sta_cnt:end_cnt),2))<mean(mean(ps_snips(2001:3000,sta_cnt:end_cnt),2))
                    [~,p(cell_cnt)]=ttest2(mean(ps_snips(1:1000,sta_cnt:end_cnt),2),mean(ps_snips(2001:3000,sta_cnt:end_cnt),2));
                else
                    p(cell_cnt)=2;
                end
                
                for jnd=1:length(mo_start)
                    try
                        if min(ROI_act_ms(mo_start(jnd)-l_win:mo_start(jnd)+r_win))>0
                            mo_snip_cnt=mo_snip_cnt+1;
                            mo_snips(:,mo_snip_cnt)=ROI_act_ms(mo_start(jnd)-l_win:mo_start(jnd)+r_win);
                            curr_fb_session=find(mo_start(jnd)>fb_onsets&mo_start(jnd)<fb_offsets);
                            if isempty(curr_fb_session)
                                curr_fb_session=0;
                            end
                            mo_snip_time(mo_snip_cnt,1)=curr_fb_session;
                            mo_snip_time(mo_snip_cnt,2)=cell_cnt;
                        end
                    end
                end
                for jnd=1:length(vs_start)
                    try
                        if min(ROI_act_ms(vs_start(jnd)-l_win:vs_start(jnd)+r_win))>0
                            vs_snip_cnt=vs_snip_cnt+1;
                            vs_snips(:,vs_snip_cnt)=ROI_act_ms(vs_start(jnd)-l_win:vs_start(jnd)+r_win);
                        end
                    end
                end
                
                for jnd=1:length(vs_stop)
                    try
                        if min(ROI_act_ms(vs_stop(jnd)-l_win:vs_stop(jnd)+r_win))>0
                            vs_stop_snip_cnt=vs_stop_snip_cnt+1;
                            vs_stop_snips(:,vs_stop_snip_cnt)=ROI_act_ms(vs_stop(jnd)-l_win:vs_stop(jnd)+r_win);
                        end
                    end
                end
                
                for jnd=1:length(saccade_times_tmp)
                    try
                        tmp=ROI_act_ms(round(saccade_times_tmp(jnd)-l_win):round(saccade_times_tmp(jnd)+r_win));
                        if min(tmp)>0
                            saccade_cnt=saccade_cnt+1;
                            saccade_snips(:,saccade_cnt)=tmp;
                         end
                    end
                end
                
                 act=AV_data(ind).raw_act(:,knd)-1;
                 act=act.*(act>(3.72*std_noise))+1;
                 %act=AV_data(ind).raw_act(:,knd);
                
%                 type_transitions=find(diff(AV_data(ind).velP_ind)==-1);
%                 
%                 for jnd=1:length(type_transitions)
%                     type(type_transitions(jnd):type_transitions(jnd)+200)=-1;
%                 end
                
                for jnd=1:4
                    act_ind(jnd,cell_cnt)=mean(act(type(round(AV_data(ind).frame_times/10))==(jnd-1)));
                    cell_id(1,cell_cnt)=ind;
                    cell_id(2,cell_cnt)=knd;
                    cell_id(3,cell_cnt)=AV_data(ind).ExpID;
                    cell_id(4,cell_cnt)=AV_data(ind).same_site_as;
                end
                
                pb_snips=[];
                use_it=[];
                for jnd=1:length(pb_onsets)
                    try
                        if max(ROI_act_ms(pb_onsets(jnd):pb_onsets(jnd)+pb_dur))>(1+5*std_noise)&sum(AV_data(ind).velP_ind(round(pb_onsets(jnd)/10):round((pb_onsets(jnd)+pb_dur)/10)))<1000
                            use_it(jnd)=1;
                        else
                            use_it(jnd)=0;
                        end
                        pb_snips(:,jnd)=ROI_act_ms(pb_onsets(jnd):pb_onsets(jnd)+pb_dur);
                    end
                end
                tmp=[];
                pb_snips=pb_snips(1:end-10000,:);
                pb_snips=double(pb_snips>(1+4.2*std_noise));
% % %                 for jnd=2:size(pb_snips,2)
% % %                     xc_tmp=xcorr(pb_snips(:,1),pb_snips(:,jnd),'biased');
% % %                     [~,tmp(jnd-1)]=max(xc_tmp(length(pb_snips)-10000:length(pb_snips)+10000));
% % %                 end
                [rccoef,pccoef]=corrcoef(pb_snips(:,1),mean(pb_snips(:,2:end)'));
                pcf(cell_cnt,1)=pccoef(2);
                pcf(cell_cnt,2)=rccoef(2);
                
% % %                 tmp=tmp-10001;
% % %                 xcorr_delay(cell_cnt)=mean(tmp(logical(use_it(2:end))));
                
% %                 if xcorr_delay(cell_cnt)<-10 & pcf(cell_cnt,2)>0 & pcf(cell_cnt,1)<0.01
% %                     keyboard
% %                 end
                
            else
                disp(['activity in cell ' num2str(AV_data(ind).ExpID) ' - ' num2str(knd) ' too large'])
            end
        end
        
        
    end
end

act_ind_dark=nan(2,size(act_ind,2));
mo_start_dark_cnt=0;
mo_stop_dark_cnt=0;
mo_start_dark_snips=zeros(l_win+r_win+1,100000);
mo_stop_dark_snips=zeros(l_win+r_win+1,100000);

for ind=1:length(AV_data)
    if AV_data(ind).layer==do_layer & AV_data(ind).ACode==10
        
        [mo_start_dark]=parse_cond_starts(AV_data(ind).velP_ind,min_dur,refract_time,1);
        [mo_stop_dark]=parse_cond_starts(~AV_data(ind).velP_ind,min_dur,refract_time,1);
        
        type=(AV_data(ind).velP_ind&AV_data(ind).velM_ind)+2*(AV_data(ind).velP_ind&~AV_data(ind).velM_ind)+3*(~AV_data(ind).velP_ind&AV_data(ind).velM_ind);
        type=sign(type);
        type_offsets=find(diff(type)==-1);
        
        for knd=1:length(type_offsets)
            type(type_offsets(knd):type_offsets(knd)+100)=1;
        end
        
        for knd=1:size(AV_data(ind).raw_act,2)
            std_noise=sqrt(sum((AV_data(ind).raw_act(AV_data(ind).raw_act(:,knd)<1,knd)-1).^2/sum(AV_data(ind).raw_act(:,knd)<1)));
            [ROI_act_ms]=map_ROI_act_to_aux_data(AV_data(ind).raw_act(:,knd),AV_data(ind).frame_times,length(AV_data(ind).PS));
            for jnd=1:length(mo_start_dark)
                try
                    if min(ROI_act_ms(mo_start_dark(jnd)-l_win:mo_start_dark(jnd)+r_win))>0
                        mo_start_dark_cnt=mo_start_dark_cnt+1;
                        mo_start_dark_snips(:,mo_start_dark_cnt)=ROI_act_ms(mo_start_dark(jnd)-l_win:mo_start_dark(jnd)+r_win);
                    end
                end
            end
            
            for jnd=1:length(mo_stop_dark)
                try
                    if min(ROI_act_ms(mo_stop_dark(jnd)-l_win:mo_stop_dark(jnd)+r_win))>0
                        mo_stop_dark_cnt=mo_stop_dark_cnt+1;
                        mo_stop_dark_snips(:,mo_stop_dark_cnt)=ROI_act_ms(mo_start_dark(jnd)-l_win:mo_start_dark(jnd)+r_win);
                    end
                end
            end
            
            
            curr_cell_id=find(cell_id(4,:)==AV_data(ind).same_site_as&cell_id(2,:)==knd);
                         act=AV_data(ind).raw_act(:,knd)-1;
                         act=act.*(act>(3.72*std_noise))+1;
            %act=AV_data(ind).raw_act(:,knd);
            for jnd=1:2
                act_ind_dark(jnd,curr_cell_id)=mean(act(type(round(AV_data(ind).frame_times/10))==(jnd-1)));
            end
            
        end
    end
end


cell_w_dark=~max(isnan(act_ind_dark));

figure;hold on;
avg_act(1,1)=nanmean(act_ind(2,cell_w_dark));
avg_act(1,2)=nanstd(act_ind(2,cell_w_dark))/sqrt(sum(~isnan(act_ind(3,cell_w_dark))));
avg_act(2,1)=nanmean(act_ind(4,cell_w_dark));
avg_act(2,2)=nanstd(act_ind(4,cell_w_dark))/sqrt(sum(~isnan(act_ind(4,cell_w_dark))));
avg_act(3,1)=nanmean(act_ind_dark(2,cell_w_dark));
avg_act(3,2)=nanstd(act_ind_dark(2,cell_w_dark))/sqrt(sum(~isnan(act_ind_dark(2,cell_w_dark))));
avg_act(4,1)=nanmean(act_ind_dark(1,cell_w_dark));
avg_act(4,2)=nanstd(act_ind_dark(1,cell_w_dark))/sqrt(sum(~isnan(act_ind_dark(1,cell_w_dark))));
avg_act(5,1)=avg_act(2,1)+avg_act(3,1)-1;
plot(avg_act(:,1),'ko','linewidth',3,'markersize',10)
errorbar(avg_act(:,1),avg_act(:,2),'k','linestyle','none','linewidth',2)


% 
% figure;plot([act_ind(2,cell_w_dark);act_ind(4,cell_w_dark);act_ind_dark(2,cell_w_dark);act_ind_dark(1,cell_w_dark)],'k')

% find(cell_id(3,:)==7440&cell_id(2,:)==40)
ac=max(act_ind)>1.01;
act_ind_ao=act_ind(:,:);

cases_list={1,2,3,4,[1 2],[1 3],[1 4],[2 3],[2 4],[3 4]};
act_type=[];
for ind=1:size(act_ind_ao,2)
    for knd=1:length(cases_list)
        act_type(knd,ind)=mean(act_ind_ao(cases_list{knd},ind))/mean(act_ind_ao(setdiff([1:4],cases_list{knd}),ind));
    end
end

[~,act_type]=max(act_type);

seg_length=[];
sort_ind=[];
for ind=1:10
    curr_inds=find(act_type==ind);
    [~,curr_order]=sort(mean(act_ind_ao(cases_list{ind},curr_inds),1));
    sort_ind(end+1:end+length(curr_inds))=curr_inds(curr_order);
    seg_length(ind)=length(curr_inds);
end

seg_length(11)=sum(~ac);
figure;
cnt=0;
for ind=1:11
    axes('position',[0.05+0.8*cnt/size(act_ind,2)+ind*0.01 0.1 0.8*seg_length(ind)/size(act_ind,2) 0.8])
    if ind<11
        imagesc((act_ind_ao(:,sort_ind(cnt+1:cnt+seg_length(ind)))));
    else
        imagesc((act_ind(:,~ac)));
    end
    set(gca,'clim',[1 1.5])
    cnt=cnt+seg_length(ind);
    axis off
end
colorbar


ps_std=std(ps_snips(:,ps_snip_time(:,4)==1)')/sqrt(snip_cnt);
ps_mean=mean(ps_snips(:,ps_snip_time(:,4)==1)');

figure;
hold on
area([0 r_win],[2 2],'facecolor',[1 1 0.7],'edgecolor',[1 1 0.7])
area([-l_win:r_win r_win:-1:-l_win],[ps_mean+ps_std flipdim(ps_mean-ps_std,2)],'facecolor',[1 1 1]*0.7,'edgecolor',[1 1 1]*0.7)
plot([-l_win:r_win],ps_mean,'k','linewidth',2)
% plot([0 0],[1.03 1.12],'k')
ylim([1.04 1.11]);
box off
xlabel('Time [ms]')
ylabel('dF/F')


mo_std=std(mo_snips(:,1:mo_snip_cnt)')/sqrt(mo_snip_cnt);
mo_mean=mean(mo_snips(:,1:mo_snip_cnt)');

figure;
hold on
area([0 r_win],[2 2],'facecolor',[1 1 1]*0.9,'edgecolor',[1 1 1]*0.9)
area([-l_win:r_win r_win:-1:-l_win],[mo_mean+mo_std flipdim(mo_mean-mo_std,2)],'facecolor',[1 1 1]*0.7,'edgecolor',[1 1 1]*0.7)
plot([-l_win:r_win],mo_mean,'k','linewidth',2)
% plot([0 0],[1.03 1.12],'k')
ylim([1.03 1.07]);
box off
xlabel('Time [ms]')
ylabel('dF/F')

vs_std=std(vs_snips(:,1:vs_snip_cnt)')/sqrt(vs_snip_cnt);
vs_mean=mean(vs_snips(:,1:vs_snip_cnt)');

vs_stop_std=std(vs_stop_snips(:,1:vs_stop_snip_cnt)')/sqrt(vs_stop_snip_cnt);
vs_stop_mean=mean(vs_stop_snips(:,1:vs_stop_snip_cnt)');

figure;
hold on
area([0 r_win],[2 2],'facecolor',[0.7 1 0.7],'edgecolor',[0.7 1 0.7])
area([-l_win:r_win r_win:-1:-l_win],[vs_mean+vs_std flipdim(vs_mean-vs_std,2)],'facecolor',[1 1 1]*0.7,'edgecolor',[1 1 1]*0.7)
plot([-l_win:r_win],vs_mean,'k','linewidth',2)
% plot([0 0],[1.03 1.12],'k')
ylim([0.995 1.01]);
box off
xlabel('Time [ms]')
ylabel('dF/F')

figure;
hold on
area([0 r_win],[2 2],'facecolor',[1 1 1]*0.9,'edgecolor',[1 1 1]*0.9)
plot([-l_win:r_win],ps_mean,'y','linewidth',2)
plot([-l_win:r_win],mo_mean,'k','linewidth',2)
plot([-l_win:r_win],vs_mean,'g','linewidth',2)
ylim([1 1.1]);
box off
xlabel('Time [ms]')
ylabel('dF/F')


figure;
hold on

area([0 r_win],[2 2],-1,'facecolor',[1 1 1]*0.9,'edgecolor',[1 1 1]*0.9)
area([-l_win:r_win r_win:-1:-l_win],[vs_stop_mean+vs_stop_std flipdim(vs_stop_mean-vs_stop_std,2)]-mean(vs_stop_mean(500:l_win)),'facecolor',[1 1 1]*0.6,'linestyle','none')
area([-l_win:r_win r_win:-1:-l_win],[ps_mean+ps_std flipdim(ps_mean-ps_std,2)]-mean(ps_mean(500:l_win)),'facecolor',[1 1 1]*0.6,'linestyle','none')
area([-l_win:r_win r_win:-1:-l_win],[mo_mean+mo_std flipdim(mo_mean-mo_std,2)]-mean(mo_mean(500:l_win)),'facecolor',[1 1 1]*0.6,'linestyle','none')
area([-l_win:r_win r_win:-1:-l_win],[vs_mean+vs_std flipdim(vs_mean-vs_std,2)]-mean(vs_mean(500:l_win)),'facecolor',[1 1 1]*0.6,'linestyle','none')
plot([-500 r_win],[0 0],'k')
plot([0 0],[-0.01 0.038],'k')
plot([-l_win:r_win],vs_stop_mean-mean(vs_stop_mean(500:l_win)),'m','linewidth',2)
plot([-l_win:r_win],ps_mean-mean(ps_mean(500:l_win)),'y','linewidth',2)
plot([-l_win:r_win],mo_mean-mean(mo_mean(500:l_win)),'k','linewidth',2)
plot([-l_win:r_win],vs_mean-mean(vs_mean(500:l_win)),'g','linewidth',2)

ylim([-0.01 0.038]);
xlim([-500 r_win])
box off
xlabel('Time [s]')
set(gca,'xtick',[-l_win:1000:r_win],'xticklabel',[-1:5])
set(gca,'ytick',[-1:1:3]*1e-2,'yticklabel',[-1:1:3])
ylabel('dF/F')




vel_snip_mean=-mean(vel_snip');
vel_snip_ste=std(vel_snip')/sqrt(vel_snip_cnt);
figure;
hold on
area([0 10000],[2 2],'facecolor',[1 1 0.7],'edgecolor',[1 1 0.7])
area([-5000 10000 10000 -5000],[mean(vel_snip_mean(1:5000))+mean(vel_snip_ste(1:5000)) mean(vel_snip_mean(1:5000))+mean(vel_snip_ste(1:5000)) mean(vel_snip_mean(1:5000))-mean(vel_snip_ste(1:5000)) mean(vel_snip_mean(1:5000))-mean(vel_snip_ste(1:5000))],'linestyle','none','facecolor',[1 1 1]*0.9)
area([-5000:10000 10000:-1:-5000],[vel_snip_mean+vel_snip_ste flipdim(vel_snip_mean-vel_snip_ste,2)],'facecolor',[1 1 1]*0.7,'edgecolor',[1 1 1]*0.7)
plot([-5000:10000],-mean(vel_snip'),'k','linewidth',2)
plot([-5000 10000],[1 1]*mean(vel_snip_mean(1:5000)),'k--','linewidth',2)
box off
ylim([0 0.025])
xlabel('Time [ms]')
ylabel('Running speed [AU]')

figure;
hold on
plot(-mean(vel_snip(:,vel_pb_session==2)')+(mean(mean(vel_snip(1:5000,vel_pb_session==2)))),'k')
plot(-mean(vel_snip(:,vel_pb_session==3)')+(mean(mean(vel_snip(1:5000,vel_pb_session==3)))),'r')
plot(-mean(vel_snip(:,vel_pb_session==4)')+(mean(mean(vel_snip(1:5000,vel_pb_session==4)))),'c')



max_ps_snips=max(ps_snips(2000:3000,1:snip_cnt))-max(ps_snips(1:1000,1:snip_cnt));
max_ps_snips_pre=mean(ps_snips(1:1000,1:snip_cnt));
avg_peak_ps_resp=[];
ste_peak_ps_resp=[];
step_size=1e5;
win_size=3e5;
for ind=1:round(1e6/step_size)-2
    l_bound=(ind-1)*step_size;
    u_bound=(ind-1)*step_size+win_size;
    avg_peak_ps_resp(ind)=mean(max_ps_snips(ps_snip_time(:,1)>l_bound&ps_snip_time(:,1)<u_bound));
    ste_peak_ps_resp(ind)=std(max_ps_snips(ps_snip_time(:,1)>l_bound&ps_snip_time(:,1)<u_bound))/sqrt(sum(ps_snip_time(:,1)>l_bound&ps_snip_time(:,1)<u_bound));
end

figure;
hold on
area([1:ind ind:-1:1],[avg_peak_ps_resp+ste_peak_ps_resp flipdim(avg_peak_ps_resp-ste_peak_ps_resp,2)],'facecolor',[1 1 1]*0.7,'linestyle','none')
plot(avg_peak_ps_resp,'k','linewidth',2)
% plot(avg_peak_ps_resp_pre,'k')
figure;
plot(ps_snip_time(:,5),max_ps_snips,'.')

figure;hold on
area([0 r_win],[2 2],-1,'facecolor',[1 1 0.7],'edgecolor',[1 1 0.7])
colors='kbrc';
for ind=3:-1:1
    tmp=(mean(ps_snips(:,ps_snip_time(:,4)==ind+1)'));
    plot([-l_win:r_win],tmp-mean(tmp(750:1000)),'color',colors(ind),'linewidth',2)
end
xlim([-500 r_win])
ylim([-0.015 0.04])
box off
xlabel('Time [ms]')
ylabel('Fluorescence change [dF/F]')



figure;hold on
area([0 r_win],[2 2],-1,'facecolor',[1 1 0.7],'edgecolor',[1 1 0.7])
for ind=1:3
    tmp=mean(ps_snips(:,ps_snip_time(:,5)<-0.01*(ind-1)+0.001&ps_snip_time(:,5)>-0.01*(ind))');
    plot([-l_win:r_win],tmp-mean(tmp(500:1000)),'color',colors(ind),'linewidth',2)
end
xlim([-500 r_win])
ylim([-0.015 0.04])
box off
xlabel('Time [ms]')
ylabel('Fluorescence change [dF/F]')






